home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 4 / MacMania 4.toast / / Demo's / Igor Demo Pro / 3 PutContentsIn Igor Pro Folder / Technical Notes / Igor Tech Notes / TN006 DSP support / Windows & PSD TEXT < prev   
Text File  |  1994-10-31  |  5KB  |  244 lines
  1.  
  2. Function log2(val)
  3.     variable val
  4.     
  5.     return ln(val)/ln(2)
  6. end
  7.  
  8.  
  9. Function CeilPwr2(x)
  10.     variable x
  11.     
  12.     return 2^(ceil(log2(x)))
  13. end
  14.  
  15.  
  16.  
  17. | brings the first graph window containing the given wave to the front.  If no such window is found
  18. | then one is created.
  19. Proc BringDestFront(w)
  20.     string w
  21.     
  22.     string win
  23.     variable winIndex,didit=0
  24.     
  25.     CheckDisplayed /A $w
  26.     if(V_flag)                    | if displayed somewhere, need to know if is graph & bring to front if so
  27.         do
  28.             win=WinName(winIndex, 1)        | name of top graph window
  29.             if( cmpstr(win,"") == 0 )
  30.                 break;                            | no more graph wndows
  31.             endif
  32.             DoWindow /F $win
  33.             CheckDisplayed  $w
  34.             if(V_Flag)
  35.                 didit= 1
  36.                 break
  37.             endif
  38.             winIndex += 1
  39.         while(1)                | exit via break
  40.     endif
  41.     if(!didit)
  42.         Display $w                | if not displayed anywhere then ok to create a new window
  43.     endif
  44. end
  45.  
  46.  
  47.  
  48.  
  49. |******** several window functions...
  50. |    each returns the average sum square value
  51. |    this is needed for psd normalization 
  52. |********
  53.  
  54.  
  55. Function ParzenOrWelch(w,isWelch)
  56.     wave w
  57.     variable isWelch
  58.     
  59.     variable N=numpnts(w),tmp
  60.  
  61.     variable sumsqr=0,jj=0
  62.     do
  63.         tmp= (jj-0.5*(N-1))/(0.5*(N+1))
  64.         if(isWelch)
  65.             tmp= 1 - tmp^2
  66.         else
  67.             tmp= 1-abs(tmp)
  68.         endif
  69.         sumsqr += tmp^2
  70.         w[jj] *= tmp
  71.         jj+=1
  72.     while(jj<N)
  73.     return sumsqr/N
  74. end
  75.  
  76. Function Parzen(w)
  77.     wave w
  78.  
  79.     return ParzenOrWelch(w,0)
  80. End
  81.  
  82. Function Welch(w)
  83.     wave w
  84.  
  85.     return ParzenOrWelch(w,1)
  86. End
  87.  
  88. Function Kaiser(w,beta)    | see D.F. Elliott, Handbook of Digital SIgnal Processing, Academic Press, P 68
  89.     wave w
  90.     variable beta
  91.     
  92.     variable N=numpnts(w),tmp
  93.     variable IBeta= BessI(0,beta), NM1= (N-1)/2
  94.  
  95.     variable sumsqr=0,jj=0
  96.     do
  97.         tmp= BessI(0,beta*sqrt(1 - ( (jj-NM1)/NM1 )^2))/IBeta
  98.         sumsqr += tmp^2
  99.         w[jj] *= tmp
  100.         jj+=1
  101.     while(jj<N)
  102.     return sumsqr/N
  103. end
  104.  
  105.  
  106. Function Hamming(w)
  107.     wave w
  108.     
  109.     variable N=numpnts(w),tmp
  110.     variable omega= 2*Pi/N,mid= (N-1)/2
  111.  
  112.     variable sumsqr=0,i=0
  113.     do
  114.         tmp= 0.54 + 0.46*cos(omega*(i-mid))
  115.         sumsqr += tmp^2
  116.         w[i] *= tmp
  117.         i+=1
  118.     while(i<N)
  119.     return sumsqr/N
  120. end
  121.  
  122. Function BlackmanHarris3(w)
  123.     wave w
  124.     
  125.     variable N=numpnts(w),tmp
  126.     variable omega= 2*Pi/N,mid= (N-1)/2
  127.  
  128.     variable sumsqr=0,i=0
  129.     do
  130.         tmp= 0.42323 + 0.49755*cos(omega*(i-mid)) +  0.07922*cos(omega*2*(i-mid))
  131.         sumsqr += tmp^2
  132.         w[i] *= tmp
  133.         i+=1
  134.     while(i<N)
  135.     return sumsqr/N
  136. end
  137.  
  138. Function KaiserBessel(w)
  139.     wave w
  140.     
  141.     variable N=numpnts(w),tmp
  142.     variable omega= 2*Pi/N,mid= (N-1)/2
  143.  
  144.     variable sumsqr=0,i=0
  145.     do
  146.         tmp= 0.40243 + 0.49804*cos(omega*(i-mid)) 
  147.         tmp += 0.09831*cos(omega*2*(i-mid)) +  0.00122*cos(omega*3*(i-mid))
  148.         sumsqr += tmp^2
  149.         w[i] *= tmp
  150.         i+=1
  151.     while(i<N)
  152.     return sumsqr/N
  153. end
  154.  
  155.  
  156.  
  157.  
  158. | Given a long data wave create a short result wave containing the Power
  159. | Spectral Density.  The name of the new wave is the name of the source + _psd
  160. | Note: if you don't want the baggage of all the window functions, just hardwire in
  161. | your favorite and remove the window parameter (or choose a subset)
  162. | Also, note that a unity amplitude sine wave will give an integral peak size of
  163. | 0.5 rather than 0.707 as one would get if a real sine wave of 1 Volt peak amplitude
  164. | were applied across a 1 ohm resistor.  Thus if you want physical meaning for
  165. | the results you should scale by 1.414/R where R is the actual resistance.
  166. |
  167. |    Change 941031: Divided DC component by 2
  168. | 
  169. Macro PSD(w,seglen,window)
  170.     string w
  171.     Prompt w "data wave:",popup WaveList("*",";","")
  172.     variable seglen=1
  173.     Prompt seglen,"segment length:",popup "256;512;1024;2048;4096;8192"
  174.     variable window=2
  175.     Prompt window,"Window type:",popup "Square;Hann;Parzen;Welch;Hamming;"
  176.                         "BlackmanHarris3;KaiserBessel"
  177. ;
  178.     PauseUpdate; silent 1
  179.     
  180.     variable npsd= 2^(7+seglen)                | number of points in group (resultant psd wave len= npsd/2+1)
  181.     variable psdOffset= npsd/2                    | offset each group by this amount
  182.     variable psdFirst=0                            | start of current group
  183.     variable nsrc= numpnts($w)
  184.     variable nsegs,winNorm                        | count of number of segements and window normalization factor
  185.     string destw=w+"_psd",srctmp=w+"_tmp"
  186.     string winw=w+"_psdWin"                    | window goes here
  187.     
  188.     if( npsd > nsrc/2 )
  189.         Abort "psd: source wave should be MUCH longer than the segment length"
  190.     endif
  191.     make/o/n=(npsd/2+1) $destw
  192.     make/o/n=(npsd) $srctmp,$winw; $winw= 1
  193.     if( window==1 )
  194.         winNorm= 1
  195.     else
  196.         if( window==2 )
  197.             Hanning $winw;winNorm=0.372                |  winNorm is avg squared value
  198.         else
  199.             if( window==3 )
  200.                 winNorm= Parzen($winw)
  201.             else
  202.                 if( window==4 )
  203.                     winNorm= Welch($winw)
  204.                 else
  205.                     if( window==5 )
  206.                         winNorm= Hamming($winw)
  207.                     else
  208.                         if( window==6 )
  209.                             winNorm= BlackmanHarris3($winw)
  210.                         else
  211.                             if( window==7 )
  212.                                 winNorm= KaiserBessel($winw)
  213.                             else
  214.                                 Abort "unknown window index"
  215.                             endif
  216.                         endif
  217.                     endif
  218.                 endif
  219.             endif
  220.         endif
  221.     endif    | (kinda makes you wish we had elseif or switch constructs, huh?)
  222.  
  223.     Duplicate/O/R=[0,npsd-1] $w $srctmp; $srctmp *= $winw; fft $srctmp
  224.     CopyScales/P $srctmp, $destw
  225.     $destw= magsqr($srctmp)
  226.     psdFirst= psdOffset
  227.     nsegs=1
  228.     do
  229.         Duplicate/O/R=[psdFirst,psdFirst+npsd-1] $w $srctmp;   $srctmp *= $winw
  230.         fft $srctmp;   $destw += magsqr($srctmp);   psdFirst += psdOffset; nsegs+=1
  231.     while( psdFirst+npsd < nsrc )
  232.     winNorm= 2/(winNorm*nsegs*npsd^2);  $destw *= winNorm
  233.     $destw[0] /= 2
  234.     
  235.     KillWaves $srctmp,$winw
  236.     BringDestFront(destw)
  237.     if( numpnts($destw) <= 129 )
  238.         Modify mode($destw)=4,marker($destw)=19,msize($destw)=1
  239.     else
  240.         Modify mode($destw)=0
  241.     endif
  242. end
  243.  
  244.